1 Plot delta from baseline and evolution of significantly changed clincal parameters

xlabels <- c("0", "3", "6", "9", "12", "15", "18", "21", "24")

full_metadata <- as_data_frame(sample_data(ps_full))

full_metadata_distinct <- full_metadata%>% # to not have multiple data entries per visit/patient 
  distinct(id_visit, .keep_all = T)%>%
  filter(!is.na(visit_cal_cor))%>%
  mutate(baseline_ppfev1 = case_when(visit_cal_9=="1"~ pp_fev_percent))%>%
  group_by(id)%>%
  fill(baseline_ppfev1)%>%
  mutate(delta_ppfev1 = pp_fev_percent-baseline_ppfev1)%>%
  ungroup()%>%
  mutate(baseline_abtt_days = case_when(visit_cal_9 == "1" ~ treatmentdays_365prior_visit))%>%
  group_by(id)%>%
  fill(baseline_abtt_days)%>%
  mutate(delta_atb_tt_days = treatmentdays_365prior_visit-baseline_abtt_days)%>%
  ungroup() %>% 
  mutate(baseline_chloride = case_when(visit_cal_9 == "1" ~ chlorid_mmol_cl_l))%>%
  group_by(id)%>%
  fill(baseline_chloride)%>%
  mutate(delta_chloride = chlorid_mmol_cl_l-baseline_chloride)%>%
  ungroup() %>% 
  mutate(baseline_igG = case_when(visit_cal_9 == "1" ~ ig_g_g_l))%>%
  group_by(id)%>%
  fill(baseline_igG)%>%
  mutate(delta_igG = ig_g_g_l-baseline_igG)%>%
  ungroup()
#### delta ppfev1 ####
lm <- summary(lmerTest::lmer(delta_ppfev1~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(delta_ppfev1)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_Dppfev <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_Dppfev $N_obs <- replace_na(lmm_Dppfev $N_obs, lmm_Dppfev $Total_N_obs[1])

lmm_Dppfev  <- lmm_Dppfev  %>% 
  mutate(Parameter = "Delta ppFEV1") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
Dppfev1<- full_metadata_distinct%>%
    ggplot(aes(y=delta_ppfev1, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id), color="black", alpha=0.5)+
    theme_classic()+
    theme_classic()+
    ylab("Delta ppFEV1")+
    xlab("Months from ETI treatment start")+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 43, step.increase = 0.05)
Dppfev1

#### delta chloride ####
lm <- summary(lmerTest::lmer(delta_chloride~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(delta_chloride)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_Dchloride <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_Dchloride $N_obs <- replace_na(lmm_Dchloride $N_obs, lmm_Dchloride $Total_N_obs[1])

lmm_Dchloride  <- lmm_Dchloride %>% 
  mutate(Parameter = "Delta sweat chloride (mmol/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
lm_stats <- lm_stats %>% 
  filter(Months_after_ETI_start%in%c("3","12","18","24"))

Dchloride<- full_metadata_distinct%>%
  filter(visit_cal_9%in%c("1","2","5","7","9")) %>% 
    ggplot(aes(y=delta_chloride, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
        geom_line(aes(group=id), color="black", alpha=0.5)+
    theme_classic()+
    theme_classic()+
    ylab("Delta sweat chloride (mmol/l)")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= c("0","3","12","18","24"))+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 4, step.increase = 0.05)
Dchloride

#### delta IgG ####
lm <- summary(lmerTest::lmer(delta_igG~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(delta_igG)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_Digg <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_Digg $N_obs <- replace_na(lmm_Digg $N_obs, lmm_Digg$Total_N_obs[1])

lmm_Digg  <- lmm_Digg %>% 
  mutate(Parameter = "Delta IgG (g/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
Digg<- full_metadata_distinct%>%
    ggplot(aes(y=delta_igG, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
        geom_line(aes(group=id), color="black", alpha=0.5)+
    theme_classic()+
    theme_classic()+
    ylab("Delta IgG (g/l)")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 2.4, step.increase = 0.05)
Digg

#### delta antibiotic treatment days ####
lm <- summary(lmerTest::lmer(delta_atb_tt_days~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(delta_atb_tt_days)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_Datb <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_Datb $N_obs <- replace_na(lmm_Datb $N_obs, lmm_Datb $Total_N_obs[1])

lmm_Datb  <- lmm_Datb%>% 
  mutate(Parameter = "Delta antibiotic tt days") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
Datb<- full_metadata_distinct%>%
    ggplot(aes(y=delta_atb_tt_days, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id), color="black", alpha=0.5)+
    theme_classic()+
    theme_classic()+
    ylab("Delta antibiotic tt days")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 23, step.increase = 0.05)
Datb

ggarrange(Dchloride, Dppfev1, Digg, Datb, nrow = 1)

lm_stats_1 <- rbind(lmm_Dppfev, lmm_Dchloride, lmm_Digg, lmm_Datb)

kable(lm_stats_1)
Parameter Months_after_ETI_start N_obs Total_N_obs Estimate Std..Error df t.value p fdr fdr_star
Delta ppFEV1 Baseline (Intercept) 35 235 1.0158422 3.4045717 38.64296 0.2983759 0.76701 0.84372 ns
Delta ppFEV1 3 30 235 12.4628627 1.5482261 193.34212 8.0497692 0.00000 0.00000 ***
Delta ppFEV1 6 31 235 11.3119118 1.5384045 194.31464 7.3530151 0.00000 0.00000 ***
Delta ppFEV1 9 31 235 13.4351520 1.5389186 194.47907 8.7302550 0.00000 0.00000 ***
Delta ppFEV1 12 26 235 9.3934309 1.6321052 196.03302 5.7554077 0.00000 0.00000 ***
Delta ppFEV1 15 22 235 10.2279562 1.7298648 196.79571 5.9125754 0.00000 0.00000 ***
Delta ppFEV1 18 20 235 7.4812868 1.7930864 198.43422 4.1722956 0.00005 0.00008 ***
Delta ppFEV1 21 22 235 6.8596951 1.7354919 198.80899 3.9525942 0.00011 0.00017 ***
Delta ppFEV1 24 18 235 6.3404820 1.8748622 200.42886 3.3818389 0.00087 0.00119 **
Delta ppFEV1 sex 235 NA -1.4902106 2.9866734 32.30319 -0.4989533 0.62119 0.75924 ns
Delta ppFEV1 age_y 235 NA -0.0100933 0.1079618 33.07653 -0.0934894 0.92608 0.92608 ns
Delta sweat chloride (mmol/l) Baseline (Intercept) 34 116 -0.4461111 4.0390632 40.67402 -0.1104492 0.91260 0.91260 ns
Delta sweat chloride (mmol/l) 3 24 116 -36.3784810 2.7750387 78.98507 -13.1091795 0.00000 0.00000 ***
Delta sweat chloride (mmol/l) 6 6 116 -34.7004282 4.9665859 90.79944 -6.9867771 0.00000 0.00000 ***
Delta sweat chloride (mmol/l) 9 3 116 -31.8304483 6.8780454 90.56285 -4.6278334 0.00001 0.00002 ***
Delta sweat chloride (mmol/l) 12 18 116 -40.9264249 3.0860935 81.88236 -13.2615635 0.00000 0.00000 ***
Delta sweat chloride (mmol/l) 15 5 116 -38.3445596 5.3359683 88.97815 -7.1860546 0.00000 0.00000 ***
Delta sweat chloride (mmol/l) 18 14 116 -44.8826099 3.4289466 83.23005 -13.0893287 0.00000 0.00000 ***
Delta sweat chloride (mmol/l) 21 2 116 -31.6546927 8.1648932 87.76237 -3.8769267 0.00020 0.00028 ***
Delta sweat chloride (mmol/l) 24 10 116 -42.4302600 3.9106492 82.87276 -10.8499274 0.00000 0.00000 ***
Delta sweat chloride (mmol/l) sex 116 NA 3.2813715 3.4512967 30.08667 0.9507648 0.34930 0.42692 ns
Delta sweat chloride (mmol/l) age_y 116 NA -0.0477268 0.1255767 31.08161 -0.3800613 0.70648 0.77713 ns
Delta IgG (g/l) Baseline (Intercept) 30 200 1.0007796 0.7210072 32.36255 1.3880300 0.17461 0.19207 ns
Delta IgG (g/l) 3 26 200 -1.2117237 0.3043869 163.41537 -3.9808664 0.00010 0.00016 ***
Delta IgG (g/l) 6 26 200 -1.5615271 0.3056767 164.34128 -5.1084264 0.00000 0.00000 ***
Delta IgG (g/l) 9 26 200 -1.6885322 0.3060615 164.68197 -5.5169710 0.00000 0.00000 ***
Delta IgG (g/l) 12 21 200 -1.7800419 0.3288477 165.79080 -5.4129674 0.00000 0.00000 ***
Delta IgG (g/l) 15 18 200 -1.7938293 0.3478857 166.52434 -5.1563747 0.00000 0.00000 ***
Delta IgG (g/l) 18 18 200 -1.5590856 0.3487345 167.71502 -4.4706945 0.00001 0.00003 ***
Delta IgG (g/l) 21 18 200 -1.6565225 0.3486516 168.24383 -4.7512256 0.00000 0.00001 ***
Delta IgG (g/l) 24 17 200 -1.2153766 0.3592677 170.20725 -3.3829275 0.00089 0.00122 **
Delta IgG (g/l) sex 200 NA 0.3700112 0.6559188 27.74904 0.5641113 0.57721 0.57721 ns
Delta IgG (g/l) age_y 200 NA -0.0498634 0.0234825 28.44814 -2.1234293 0.04254 0.05200 .
Delta antibiotic tt days Baseline (Intercept) 35 236 12.6828934 12.6916885 44.99841 0.9993070 0.32299 0.32299 ns
Delta antibiotic tt days 3 29 236 -11.3029382 7.3110712 197.07435 -1.5460030 0.12371 0.13608 ns
Delta antibiotic tt days 6 30 236 -16.4788730 7.2557827 197.94227 -2.2711365 0.02421 0.03330 *
Delta antibiotic tt days 9 31 236 -29.4923614 7.1817520 197.64750 -4.1065692 0.00006 0.00011 ***
Delta antibiotic tt days 12 25 236 -43.9172016 7.7131424 199.65497 -5.6938145 0.00000 0.00000 ***
Delta antibiotic tt days 15 23 236 -42.1089966 7.9473336 200.16922 -5.2985062 0.00000 0.00000 ***
Delta antibiotic tt days 18 22 236 -46.5567583 8.0877718 201.12244 -5.7564382 0.00000 0.00000 ***
Delta antibiotic tt days 21 22 236 -43.8305236 8.0833203 201.42531 -5.4223416 0.00000 0.00000 ***
Delta antibiotic tt days 24 19 236 -47.3194488 8.5588867 202.58688 -5.5286921 0.00000 0.00000 ***
Delta antibiotic tt days sex 236 NA 23.4806043 10.8309936 34.14058 2.1679086 0.03723 0.04550 *
Delta antibiotic tt days age_y 236 NA -1.0017944 0.3929336 34.92744 -2.5495261 0.01533 0.02410 *
#write_csv(lm_stats_1, "/Users/rebecca/Documents/Forschung/IMMProveCF/R_analysis/table_outputs/lmm_delta_clinical_par.csv")

2 combine plots for Figure 1

phylum_bars <- readRDS("/Users/rebecca/Documents/Forschung/IMMProveCF/R_analysis/figures/samples_phylum_wControls.rds")

library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
# Arrange the ggplot objects in a grid
grid_arranged <- grid.arrange(
  phylum_bars + labs(tag = "A") + theme(legend.margin = margin(0, 0, 0, 0), legend.box.spacing = unit(0, units = "pt"), legend.box.margin = margin(0, 0, 0, 0)),
  Dchloride + labs(tag = "B"),
  Dppfev1 + labs(tag = "C"),
  Digg + labs(tag = "D"),
  Datb+ labs(tag = "E"),
  ncol = 4, nrow = 2,
  heights = c(1.5, 1),
  widths=c(1.15,1.5,1.5,1.5),
  layout_matrix = rbind(c(1,1,1,1),
                        c(2,3,4,5))
)

# Print or save the arranged plot
grid_arranged
## TableGrob (2 x 4) "arrange": 5 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-4) arrange gtable[layout]
## 2 2 (2-2,1-1) arrange gtable[layout]
## 3 3 (2-2,2-2) arrange gtable[layout]
## 4 4 (2-2,3-3) arrange gtable[layout]
## 5 5 (2-2,4-4) arrange gtable[layout]
#ggsave("/Users/rebecca/Documents/Forschung/IMMProveCF/Manuscript_16S/Submission_CHM/figure1.tiff",grid_arranged, dpi = 600, width = 17, height = 10)

3 Plot evolution of clinical parameters and run linear mixed effect models for statistical testing

#### Sweat chloride ####
lm <- summary(lmerTest::lmer(chlorid_mmol_cl_l~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(chlorid_mmol_cl_l)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_chloride <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_chloride $N_obs <- replace_na(lmm_chloride $N_obs, lmm_chloride $Total_N_obs[1])

lmm_chloride  <- lmm_chloride %>% 
  mutate(Parameter = "Sweat chloride (mmol/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
Chloride<- full_metadata_distinct%>%
    ggplot(aes(y=chlorid_mmol_cl_l, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.5))+
    theme_classic()+
    theme_classic()+
    ylab("Sweat chloride (mmol/l)")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 90, step.increase = 0.05)
Chloride

### ppfev1 ###

lm <- summary(lmerTest::lmer(pp_fev_percent~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(pp_fev_percent)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_ppfev <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_ppfev $N_obs <- replace_na(lmm_ppfev $N_obs, lmm_ppfev $Total_N_obs[1])

lmm_ppfev  <- lmm_ppfev  %>% 
  mutate(Parameter = "ppFEV1") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
ppfev1<- full_metadata_distinct%>%
    ggplot(aes(y=pp_fev_percent, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("ppFEV1")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 145, step.increase = 0.05)
ppfev1

### ppfev1 ###

lm <- summary(lmerTest::lmer(pp_fvc_percent~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(pp_fvc_percent)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_ppfvc  <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_ppfvc$N_obs <- replace_na(lmm_ppfvc$N_obs, lmm_ppfvc$Total_N_obs[1])

lmm_ppfvc  <- lmm_ppfvc   %>% 
  mutate(Parameter = "ppFEV1") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
ppfvc<- full_metadata_distinct%>%
    ggplot(aes(y=pp_fvc_percent, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("ppFVC")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 145, step.increase = 0.05)
ppfvc

#### IgG ####
lm <- summary(lmerTest::lmer(ig_g_g_l~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(ig_g_g_l)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_igg <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_igg$N_obs <- replace_na(lmm_igg$N_obs, lmm_igg$Total_N_obs[1])

lmm_igg <- lmm_igg%>% 
  mutate(Parameter = "IgG (g/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
igG <- full_metadata_distinct%>%
    ggplot(aes(y=ig_g_g_l, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("IgG (g/l)")+
    xlab("Months from ETI treatment start")+
    #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,100,200))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 25, step.increase = 0.05)
igG

#### Interleukin 6 ####

lm <- summary(lmerTest::lmer(interleukin6_pg_ml~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(interleukin6_pg_ml)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_il6 <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_il6$N_obs <- replace_na(lmm_il6$N_obs, lmm_il6$Total_N_obs[1])

lmm_il6  <- lmm_il6  %>% 
  mutate(Parameter = "Interleukin 6 (pg/ml)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
Il6 <- full_metadata_distinct%>%
    ggplot(aes(y=interleukin6_pg_ml, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("Interleukin 6 (pg/ml)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 2.1, step.increase = 0.05)
Il6

#### Interleukin 8 ####

lm <- summary(lmerTest::lmer(interleukin8_pg_ml~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(interleukin8_pg_ml)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_il8 <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_il8$N_obs <- replace_na(lmm_il8$N_obs, lmm_il8$Total_N_obs[1])

lmm_il8  <- lmm_il8  %>% 
  mutate(Parameter = "Interleukin 8 (pg/ml)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
Il8 <- full_metadata_distinct%>%
    ggplot(aes(y=interleukin8_pg_ml, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("Interleukin 8 (pg/ml)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 2.1, step.increase = 0.05)
Il8

#### CRP ####

lm <- summary(lmerTest::lmer(crp_mg_l~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(crp_mg_l)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_crp <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_crp$N_obs <- replace_na(lmm_crp$N_obs, lmm_crp$Total_N_obs[1])

lmm_crp  <- lmm_crp  %>% 
  mutate(Parameter = "CRP (mg/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
CRP <- full_metadata_distinct%>%
    ggplot(aes(y=crp_mg_l, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("CRP (mg/l)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 2.1, step.increase = 0.05)
CRP

#### Calprotectin in stool ####

lm <- summary(lmerTest::lmer(calprotectin_amg_g_stuhl~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(calprotectin_amg_g_stuhl)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_calp <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_calp$N_obs <- replace_na(lmm_calp$N_obs, lmm_calp$Total_N_obs[1])

lmm_calp <- lmm_calp %>% 
  mutate(Parameter = "Calprotectin (mg/g)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
calp_stool <- full_metadata_distinct%>%
    ggplot(aes(y=calprotectin_amg_g_stuhl, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("Calprotectin in stool (mg/g)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 3, step.increase = 0.05)
calp_stool

#### GOT/ASAT liver enzyme ####
lm <- summary(lmerTest::lmer(got_asat_u_l~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(got_asat_u_l)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_got <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_got$N_obs <- replace_na(lmm_got$N_obs, lmm_got$Total_N_obs[1])

lmm_got <- lmm_got %>% 
  mutate(Parameter = "GOT/ASAT (U/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
got<- full_metadata_distinct%>%
    ggplot(aes(y=got_asat_u_l, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("GOT/ASAT (ul/l)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 2, step.increase = 0.05)
got

#### GPT/ALAT liver enzyme ####

lm <- summary(lmerTest::lmer(gpt_alat_u_l~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(gpt_alat_u_l)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_gpt <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_gpt$N_obs <- replace_na(lmm_gpt$N_obs, lmm_gpt$Total_N_obs[1])

lmm_gpt <- lmm_gpt %>% 
  mutate(Parameter = "GPT/ALAT (U/l)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
gpt<- full_metadata_distinct%>%
    ggplot(aes(y=gpt_alat_u_l, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("GPT/ALAT (ul/l)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 2, step.increase = 0.025)
gpt

### HbA1c ###

lm <- summary(lmerTest::lmer(hb_a1c_percent~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(hb_a1c_percent)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24")) %>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_hba1c <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_hba1c$N_obs <- replace_na(lmm_hba1c$N_obs, lmm_hba1c$Total_N_obs[1])

lmm_hba1c <- lmm_hba1c %>% 
  mutate(Parameter = "HbA1c (%)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
hba1c<- full_metadata_distinct%>%
    ggplot(aes(y=hb_a1c_percent, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("HbA1c (%)")+
    xlab("Months from ETI treatment start")+
    scale_y_log10()+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 1.1, step.increase = 0.025)
hba1c

### Leucocytes ###

lm <- summary(lmerTest::lmer(leukozyten_nl~ visit_cal_9 + sex + age_y+ (1|id), data=full_metadata_distinct))

  coefs <- data.frame(coef(lm))
  fdr <- p.adjust(coefs$Pr...t.., method = "fdr", n = nrow(coefs))
  
  lm_stats <- bind_cols(coefs, fdr)%>%
    mutate(p = Pr...t..)%>%
    mutate(fdr = ...6)%>%
    select(-c(Pr...t.., ...6))%>%
    rownames_to_column()%>%
    mutate(Months_after_ETI_start=rowname)%>%
    mutate(Months_after_ETI_start= case_when(Months_after_ETI_start=="(Intercept)"~"Baseline (Intercept)", Months_after_ETI_start=="visit_cal_92"~"3", Months_after_ETI_start=="visit_cal_93"~"6",Months_after_ETI_start=="visit_cal_94"~"9",Months_after_ETI_start=="visit_cal_95"~"12",Months_after_ETI_start=="visit_cal_96"~"15",Months_after_ETI_start=="visit_cal_97"~"18",Months_after_ETI_start=="visit_cal_98"~"21",Months_after_ETI_start=="visit_cal_99"~"24", Months_after_ETI_start=="sex2"~"sex", Months_after_ETI_start=="age_y"~"age_y"))%>%
    mutate(fdr_star = case_when(fdr<=0.001 ~ "***", fdr<=0.01 ~ "**", fdr<=0.05 ~ "*",fdr<=0.1 ~ ".",fdr>=0.1 ~ "ns"))%>%
    select(-rowname)%>%
    select(Months_after_ETI_start, Estimate, Std..Error, df, t.value, p, fdr, fdr_star)%>%
    mutate(p=round(p,5))%>%
    mutate(fdr=round(fdr, 5))
## New names:
## • `` -> `...6`
lm_stats
# print lm stats into data frame and add N of observations

# print lm stats into data frame and add N of observations
nobs_per_visit <- full_metadata_distinct %>%
  filter(!is.na(leukozyten_nl)) %>% 
  group_by(visit_cal_9) %>%
  summarise(N_obs = n()) %>% 
  mutate(Months_after_ETI_start= case_when(visit_cal_9=="1"~"Baseline (Intercept)", visit_cal_9=="2"~"3", visit_cal_9=="3"~"6",visit_cal_9=="4"~"9",visit_cal_9=="5"~"12",visit_cal_9=="6"~"15",visit_cal_9=="7"~"18",visit_cal_9=="8"~"21",visit_cal_9=="9"~"24"))%>% 
  mutate(Total_N_obs = sum(N_obs)) 

# Merge the two data frames
lmm_leukos <- left_join(lm_stats, nobs_per_visit, by = "Months_after_ETI_start") 
# Fill NAs in N_obs with the resulting rowsum of N_obs
lmm_leukos$N_obs <- replace_na(lmm_leukos$N_obs, lmm_leukos$Total_N_obs[1])

lmm_leukos <- lmm_leukos %>% 
  mutate(Parameter = "Leukocytes (/nl)") %>% 
  select(Parameter, Months_after_ETI_start, N_obs, Total_N_obs, everything()) %>% 
  select(-visit_cal_9)

# to make the significant asterixes printable on the plot I have to adjust the table a little
group1 <- c("(Intercept)","1","1","1","1","1","1","1","1","1","1")
  
lm_stats <-  lm_stats%>%
    bind_cols(group1)%>%
    mutate(group1=...9)%>%
  mutate(group2=case_when(Months_after_ETI_start=="3"~ "2", Months_after_ETI_start=="6"~ "3", Months_after_ETI_start=="9"~ "4",Months_after_ETI_start=="12"~ "5", Months_after_ETI_start=="15"~ "6", Months_after_ETI_start=="18"~ "7",Months_after_ETI_start=="21"~ "8", Months_after_ETI_start=="24"~ "9"))%>%
  filter(group1!="(Intercept)") %>% 
  filter(!is.na(group2))
## New names:
## • `` -> `...9`
leucos<- full_metadata_distinct%>%
    ggplot(aes(y=leukozyten_nl, x=visit_cal_9))+
     geom_boxplot(aes(fill=visit_cal_9), outlier.shape = NA, alpha=0.7)+
    geom_point(color="black", alpha=0.7, width = 0.1)+
    scale_fill_manual(values=visit_palette)+
    geom_line(aes(group=id, alpha=0.7))+
    theme_classic()+
    theme_classic()+
    ylab("Leucocytes (/nl)")+
    xlab("Months from ETI treatment start")+
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(1,10,50,100,300))+
    scale_x_discrete(labels= xlabels)+
    theme(text=element_text(size = 14), legend.position = "none")+
    #guides(fill='none') +
   stat_pvalue_manual(lm_stats, label = "fdr_star",y.position = 1.4, step.increase = 0.05)
leucos

ggarrange(ppfev1, ppfvc, Chloride, igG, Il6, calp_stool, gpt, hba1c, nrow = 2, ncol = 4)

#ggsave("/Users/rebecca/Documents/Forschung/IMMProveCF/Manuscript_16S/Submission_CHM/supp_figure1.tiff", width = 14, height = 7)
stats_table <- rbind(lmm_chloride,lmm_Dppfev,lmm_ppfvc,lmm_igg,lmm_leukos, lmm_il6, lmm_il8, lmm_calp, lmm_gpt, lmm_got, lmm_hba1c)
#write_csv(stats_table, "/Users/rebecca/Documents/Forschung/IMMProveCF/R_analysis/table_outputs/lmm_clinical_parameters.csv")
kable(stats_table)
Parameter Months_after_ETI_start N_obs Total_N_obs Estimate Std..Error df t.value p fdr fdr_star
Sweat chloride (mmol/l) Baseline (Intercept) 34 119 69.0854861 4.3279602 42.47775 15.9625974 0.00000 0.00000 ***
Sweat chloride (mmol/l) 3 25 119 -35.6683953 2.7298654 82.09485 -13.0659903 0.00000 0.00000 ***
Sweat chloride (mmol/l) 6 6 119 -35.0235154 4.9688991 91.98396 -7.0485463 0.00000 0.00000 ***
Sweat chloride (mmol/l) 9 3 119 -37.3910610 6.8785114 91.54992 -5.4359234 0.00000 0.00000 ***
Sweat chloride (mmol/l) 12 19 119 -40.6820780 3.0234061 84.78035 -13.4557106 0.00000 0.00000 ***
Sweat chloride (mmol/l) 15 5 119 -39.2537245 5.3309544 90.11204 -7.3633578 0.00000 0.00000 ***
Sweat chloride (mmol/l) 18 14 119 -44.5773069 3.4142706 85.44003 -13.0561726 0.00000 0.00000 ***
Sweat chloride (mmol/l) 21 2 119 -32.9562276 8.1513534 88.86314 -4.0430375 0.00011 0.00014 ***
Sweat chloride (mmol/l) 24 11 119 -42.0315021 3.7597756 86.30286 -11.1792582 0.00000 0.00000 ***
Sweat chloride (mmol/l) sex 119 NA 10.0339457 3.6951540 32.75024 2.7154337 0.01048 0.01048 *
Sweat chloride (mmol/l) age_y 119 NA 0.4255756 0.1348529 33.68501 3.1558499 0.00336 0.00370 **
Delta ppFEV1 Baseline (Intercept) 35 235 1.0158422 3.4045717 38.64296 0.2983759 0.76701 0.84372 ns
Delta ppFEV1 3 30 235 12.4628627 1.5482261 193.34212 8.0497692 0.00000 0.00000 ***
Delta ppFEV1 6 31 235 11.3119118 1.5384045 194.31464 7.3530151 0.00000 0.00000 ***
Delta ppFEV1 9 31 235 13.4351520 1.5389186 194.47907 8.7302550 0.00000 0.00000 ***
Delta ppFEV1 12 26 235 9.3934309 1.6321052 196.03302 5.7554077 0.00000 0.00000 ***
Delta ppFEV1 15 22 235 10.2279562 1.7298648 196.79571 5.9125754 0.00000 0.00000 ***
Delta ppFEV1 18 20 235 7.4812868 1.7930864 198.43422 4.1722956 0.00005 0.00008 ***
Delta ppFEV1 21 22 235 6.8596951 1.7354919 198.80899 3.9525942 0.00011 0.00017 ***
Delta ppFEV1 24 18 235 6.3404820 1.8748622 200.42886 3.3818389 0.00087 0.00119 **
Delta ppFEV1 sex 235 NA -1.4902106 2.9866734 32.30319 -0.4989533 0.62119 0.75924 ns
Delta ppFEV1 age_y 235 NA -0.0100933 0.1079618 33.07653 -0.0934894 0.92608 0.92608 ns
ppFEV1 Baseline (Intercept) 35 234 101.1293680 5.8076667 34.29655 17.4130806 0.00000 0.00000 ***
ppFEV1 3 30 234 10.6705013 1.5983305 190.94312 6.6760293 0.00000 0.00000 ***
ppFEV1 6 31 234 9.5525085 1.5907892 192.46793 6.0048864 0.00000 0.00000 ***
ppFEV1 9 31 234 11.6134594 1.5922728 193.23643 7.2936368 0.00000 0.00000 ***
ppFEV1 12 26 234 9.4728823 1.6929252 195.44954 5.5955706 0.00000 0.00000 ***
ppFEV1 15 21 234 8.2862901 1.8252527 197.22898 4.5398043 0.00001 0.00002 ***
ppFEV1 18 20 234 6.7752830 1.8704931 201.28564 3.6221908 0.00037 0.00051 ***
ppFEV1 21 22 234 6.7483801 1.8122199 202.30180 3.7238197 0.00025 0.00040 ***
ppFEV1 24 18 234 4.9188976 1.9654334 205.93999 2.5027038 0.01310 0.01601 *
ppFEV1 sex 234 NA 5.4206413 5.2596917 31.46821 1.0306006 0.31058 0.31058 ns
ppFEV1 age_y 234 NA -0.3240456 0.1884772 33.01347 -1.7192824 0.09493 0.10442 ns
IgG (g/l) Baseline (Intercept) 30 232 10.2534962 1.2180677 34.63855 8.4178378 0.00000 0.00000 ***
IgG (g/l) 3 29 232 -1.2524893 0.2837596 189.26840 -4.4139106 0.00002 0.00002 ***
IgG (g/l) 6 31 232 -1.7204948 0.2818769 191.18140 -6.1037089 0.00000 0.00000 ***
IgG (g/l) 9 31 232 -1.8475817 0.2825581 192.28914 -6.5387675 0.00000 0.00000 ***
IgG (g/l) 12 26 232 -1.9040301 0.3005890 195.01382 -6.3343315 0.00000 0.00000 ***
IgG (g/l) 15 22 232 -1.9622117 0.3183205 197.82439 -6.1642640 0.00000 0.00000 ***
IgG (g/l) 18 22 232 -1.7266318 0.3208876 202.32241 -5.3807989 0.00000 0.00000 ***
IgG (g/l) 21 22 232 -1.8791860 0.3214042 203.76324 -5.8467995 0.00000 0.00000 ***
IgG (g/l) 24 19 232 -1.5403600 0.3402355 208.91737 -4.5273346 0.00001 0.00001 ***
IgG (g/l) sex 232 NA -1.3508350 1.1103225 31.91779 -1.2166150 0.23267 0.23267 ns
IgG (g/l) age_y 232 NA 0.1321313 0.0395804 34.22812 3.3383045 0.00204 0.00225 **
Leukocytes (/nl) Baseline (Intercept) 32 234 9.3385414 1.2322962 73.25221 7.5781629 0.00000 0.00000 ***
Leukocytes (/nl) 3 30 234 -1.8968711 1.1078251 199.40245 -1.7122478 0.08841 0.20730 ns
Leukocytes (/nl) 6 30 234 -2.0856793 1.1105841 201.01311 -1.8780021 0.06183 0.20730 ns
Leukocytes (/nl) 9 30 234 -1.7918893 1.1099385 200.76492 -1.6144042 0.10801 0.20730 ns
Leukocytes (/nl) 12 26 234 -1.7316513 1.1595152 203.45800 -1.4934270 0.13687 0.20730 ns
Leukocytes (/nl) 15 23 234 -1.6622048 1.2059396 203.11072 -1.3783484 0.16961 0.20730 ns
Leukocytes (/nl) 18 22 234 0.6454324 1.2252530 204.03064 0.5267748 0.59892 0.65881 ns
Leukocytes (/nl) 21 22 234 -1.7581928 1.2239677 204.27608 -1.4364699 0.15240 0.20730 ns
Leukocytes (/nl) 24 19 234 -1.8208489 1.2918426 205.88786 -1.4094975 0.16020 0.20730 ns
Leukocytes (/nl) sex 234 NA -1.7000438 0.9172221 34.43945 -1.8534701 0.07240 0.20730 ns
Leukocytes (/nl) age_y 234 NA -0.0004959 0.0336626 36.37515 -0.0147317 0.98833 0.98833 ns
Interleukin 6 (pg/ml) Baseline (Intercept) 23 142 12.0913248 3.3496209 66.71004 3.6097591 0.00059 0.00527 **
Interleukin 6 (pg/ml) 3 17 142 -10.4236933 3.1643478 119.79613 -3.2941048 0.00130 0.00527 **
Interleukin 6 (pg/ml) 6 14 142 -8.2930616 3.3910124 123.08216 -2.4456005 0.01588 0.03283 *
Interleukin 6 (pg/ml) 9 19 142 -10.0901876 3.0916936 118.85769 -3.2636442 0.00144 0.00527 **
Interleukin 6 (pg/ml) 12 16 142 -7.7979661 3.2482240 119.60932 -2.4006861 0.01791 0.03283 *
Interleukin 6 (pg/ml) 15 7 142 -8.3825303 4.3196965 120.45984 -1.9405368 0.05465 0.06876 .
Interleukin 6 (pg/ml) 18 12 142 -6.6609438 3.5554994 120.25715 -1.8734200 0.06344 0.06978 .
Interleukin 6 (pg/ml) 21 19 142 -5.9474023 3.0852043 119.86489 -1.9277175 0.05626 0.06876 .
Interleukin 6 (pg/ml) 24 15 142 -9.9125157 3.3380392 118.40559 -2.9695624 0.00361 0.00993 **
Interleukin 6 (pg/ml) sex 142 NA -5.5781142 2.4148587 31.84047 -2.3099133 0.02753 0.04326 *
Interleukin 6 (pg/ml) age_y 142 NA 0.1478922 0.0894327 38.45241 1.6536708 0.10634 0.10634 ns
Interleukin 8 (pg/ml) Baseline (Intercept) 28 212 33.4005300 12.6078088 101.07647 2.6491939 0.00936 0.10301 ns
Interleukin 8 (pg/ml) 3 24 212 -13.6132223 13.6578198 183.33274 -0.9967347 0.32021 0.62451 ns
Interleukin 8 (pg/ml) 6 29 212 -13.4123485 13.0144032 183.94510 -1.0305773 0.30409 0.62451 ns
Interleukin 8 (pg/ml) 9 28 212 -10.1635654 13.1531518 185.83268 -0.7727095 0.44068 0.62451 ns
Interleukin 8 (pg/ml) 12 23 212 -0.8228657 13.9117605 188.35988 -0.0591489 0.95290 0.95290 ns
Interleukin 8 (pg/ml) 15 21 212 12.5951739 14.3039541 188.67320 0.8805379 0.37969 0.62451 ns
Interleukin 8 (pg/ml) 18 21 212 -1.2850447 14.3494815 189.41538 -0.0895534 0.92874 0.95290 ns
Interleukin 8 (pg/ml) 21 22 212 10.5982778 14.1311375 189.24140 0.7499947 0.45419 0.62451 ns
Interleukin 8 (pg/ml) 24 16 212 -12.9644453 15.6343884 191.24249 -0.8292263 0.40801 0.62451 ns
Interleukin 8 (pg/ml) sex 212 NA 7.4896144 8.3854131 33.85495 0.8931718 0.37807 0.62451 ns
Interleukin 8 (pg/ml) age_y 212 NA -0.1625423 0.3145589 37.69665 -0.5167308 0.60836 0.74355 ns
Calprotectin (mg/g) Baseline (Intercept) 18 120 178.7703785 25.9977966 50.57516 6.8763665 0.00000 0.00000 ***
Calprotectin (mg/g) 3 9 120 -117.3008502 34.9875511 104.56807 -3.3526453 0.00112 0.00153 **
Calprotectin (mg/g) 6 14 120 -106.2924317 30.2901786 97.20218 -3.5091385 0.00068 0.00107 **
Calprotectin (mg/g) 9 15 120 -88.1664265 30.0300230 92.32422 -2.9359427 0.00420 0.00513 **
Calprotectin (mg/g) 12 11 120 -126.6824229 32.8491393 97.59591 -3.8564914 0.00021 0.00038 ***
Calprotectin (mg/g) 15 11 120 -139.7518941 33.0813239 97.17936 -4.2244952 0.00005 0.00012 ***
Calprotectin (mg/g) 18 12 120 -143.6487962 32.1672533 105.83157 -4.4656842 0.00002 0.00006 ***
Calprotectin (mg/g) 21 16 120 -159.0050397 29.5171626 99.20808 -5.3868674 0.00000 0.00000 ***
Calprotectin (mg/g) 24 14 120 -164.8620305 30.9167393 100.48976 -5.3324521 0.00000 0.00000 ***
Calprotectin (mg/g) sex 120 NA -12.0107707 18.4637036 19.30519 -0.6505071 0.52303 0.52303 ns
Calprotectin (mg/g) age_y 120 NA 0.9174243 0.6591936 17.76059 1.3917374 0.18119 0.19931 ns
GPT/ALAT (U/l) Baseline (Intercept) 33 236 12.8189870 5.6236280 44.68613 2.2794870 0.02747 0.04316 *
GPT/ALAT (U/l) 3 30 236 11.6287612 3.4478332 195.46765 3.3727737 0.00090 0.00445 **
GPT/ALAT (U/l) 6 31 236 10.7105478 3.4252042 196.50848 3.1269809 0.00203 0.00559 **
GPT/ALAT (U/l) 9 31 236 7.9681772 3.4261340 196.56343 2.3257051 0.02105 0.03860 *
GPT/ALAT (U/l) 12 26 236 11.9048006 3.6264203 198.53358 3.2827967 0.00121 0.00445 **
GPT/ALAT (U/l) 15 23 236 9.4939656 3.7839403 199.08265 2.5090157 0.01290 0.02839 *
GPT/ALAT (U/l) 18 21 236 6.4525604 3.9080376 200.04581 1.6510999 0.10029 0.13789 ns
GPT/ALAT (U/l) 21 22 236 3.2229866 3.8473249 200.33394 0.8377215 0.40319 0.40319 ns
GPT/ALAT (U/l) 24 19 236 6.2267434 4.0710685 201.54846 1.5295108 0.12771 0.15608 ns
GPT/ALAT (U/l) sex 236 NA 18.2526070 4.7434809 32.26636 3.8479352 0.00053 0.00445 **
GPT/ALAT (U/l) age_y 236 NA 0.1657523 0.1723532 33.22639 0.9617012 0.34315 0.37746 ns
GOT/ASAT (U/l) Baseline (Intercept) 33 236 27.0033788 2.9279484 42.07908 9.2226281 0.00000 0.00000 ***
GOT/ASAT (U/l) 3 30 236 3.1357365 1.6487210 194.82588 1.9019207 0.05866 0.21508 ns
GOT/ASAT (U/l) 6 31 236 2.7475023 1.6385718 195.89342 1.6767665 0.09518 0.26175 ns
GOT/ASAT (U/l) 9 31 236 1.6284189 1.6390431 195.97009 0.9935180 0.32168 0.58975 ns
GOT/ASAT (U/l) 12 26 236 1.7633046 1.7357418 197.80044 1.0158795 0.31093 0.58975 ns
GOT/ASAT (U/l) 15 22 236 1.4495953 1.8381258 198.49184 0.7886268 0.43127 0.59300 ns
GOT/ASAT (U/l) 18 22 236 0.5760887 1.8435197 199.47163 0.3124939 0.75499 0.86463 ns
GOT/ASAT (U/l) 21 22 236 0.5008440 1.8424279 199.80357 0.2718391 0.78603 0.86463 ns
GOT/ASAT (U/l) 24 19 236 -0.0988094 1.9502776 201.11846 -0.0506643 0.95964 0.95964 ns
GOT/ASAT (U/l) sex 236 NA 9.7185740 2.5047824 32.00660 3.8800072 0.00049 0.00269 **
GOT/ASAT (U/l) age_y 236 NA -0.0800732 0.0909004 32.90946 -0.8808894 0.38477 0.59300 ns
HbA1c (%) Baseline (Intercept) 16 123 5.9435974 0.5022168 39.10772 11.8347242 0.00000 0.00000 ***
HbA1c (%) 3 19 123 -0.6201108 0.1799150 83.35358 -3.4466875 0.00089 0.00490 **
HbA1c (%) 6 12 123 -0.6172339 0.2286378 87.89141 -2.6996145 0.00833 0.03053 *
HbA1c (%) 9 12 123 -0.5039793 0.2060707 83.87834 -2.4456628 0.01655 0.03641 *
HbA1c (%) 12 11 123 -0.5162835 0.2092032 83.60846 -2.4678569 0.01563 0.03641 *
HbA1c (%) 15 9 123 -0.4156051 0.2268026 83.73601 -1.8324528 0.07044 0.11069 ns
HbA1c (%) 18 10 123 -0.4834953 0.2131365 84.69849 -2.2684777 0.02584 0.04738 *
HbA1c (%) 21 19 123 -0.2903652 0.1779156 85.67842 -1.6320387 0.10634 0.14622 ns
HbA1c (%) 24 15 123 -0.0533455 0.1862931 86.86123 -0.2863524 0.77529 0.77529 ns
HbA1c (%) sex 123 NA -0.3215008 0.4365231 32.20279 -0.7365035 0.46676 0.51343 ns
HbA1c (%) age_y 123 NA 0.0231887 0.0158975 33.78674 1.4586412 0.15390 0.18810 ns